SUBROUTINE mayorsolution(nn,xx,fval)
  USE commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'
  
  INTEGER, INTENT(IN)  :: nn
  REAL(8), INTENT(IN)  :: xx(nn)
  REAL(8), INTENT(OUT) :: fval
  REAL(8)              :: zzprl, fundsl, popsizel, totresl, budget, util, vemax,pcons,zzpu,qcons,ln_relqcons, newsav, newsteal, new_fr_stolen, &
                          emaxp, emaxp1, emaxp2(Nfines), emaxp3(Nfines), qcons_w_riv,&
                          netW(Nfines), fine_fun, ipaul, v_slope, multi_slope, betat1_t(MMg), betat2_t(MMg), betat3_t(Nwealth)
  INTEGER              :: imunl, iedul, iterml, ub_sav, ub_netW(Nfines), isav, timel, ifine

  EXTERNAL prodfunc, utilfun_mayor
  
  !Internalizing global variables
  imunl = imunt
  iedul = iedut
  zzprl = zzprt
  popsizel = popsizeg
  fundsl = fundst
  totresl = totrest
  ipaul = ipaug
  iterml = itermg
  timel = timeg

  !Initialize variables
  util = -infty
  vemax = -infty
  pcons = -infty
  zzpu = -infty
  qcons = -infty
  ln_relqcons = -infty
  emaxp = -infty

  !Transform xx into relevant variables
  IF(flag_sopt_sa == 0) THEN
     IF (xx(1) < 10d0) THEN
        newsteal = (EXP(xx(1))/(1d0 + EXP(xx(1))))*fundsl
     ELSE
        newsteal = fundsl-10d0*small_no
     END IF
     budget = totresl + newsteal
     newsav = (budget*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2)))
  ELSE
     newsteal = xx(1)
     newsav = xx(2)
  END IF

  !Calculating utility
  pcons = totresl - newsav + newsteal

  IF (pcons <= 0.0d0) THEN
!!$     WRITE(*,*) 'ERROR in mayorsolution pcons/qcons',pcons,qcons
     util = bignumber
     fval = infty
!!$     WRITE(*,*) 'fval-inf',fval
     RETURN
  END IF

  !We use per-capita public inputs                                                                       |
  zzpu = (fundsl - newsteal)/popsizeg
  IF (zzpu <= small_no) THEN
     qcons = subqcons
     GOTO 671
  END IF

  !The output qcons is in per-capita and we account for the degree of rivalry in the utility function
  CALL prodfunc(imunl,iedul,zzpu,zzprl,qcons)

  !We transform per-capita public consumption in total public consumption
  !Below we take into account the degree of rivalry by dividing by popsize**eta
671 qcons = qcons*popsizeg

  !We are making sure that public consumption is not negative; In the simulations at the optimum it is never negative
  IF (qcons <= 0d0) qcons = small_no
  
  !We transform the level of public consumption to take into account the degree of rivalry
  !This is not the public consumption we output
  qcons_w_riv = qcons/(popsizeg**eta)

  CALL utilfun_mayor(pcons,qcons_w_riv,theta,util)

  !Calculating the future value function depends on the history of mayors -- really only relevant for iterm =2 (those who have a history)
  !Case 1: iterm = 1 & iterm = 2 with ipast = 1 (not audited mayors): Construct value function with the three different betas
  !Case 2: iterm = 2 & ipast = 2 (audited that did not steal): Construct value function only with beta2 and beta3 --they may be audited and caught stealing in the future
  !        if they don't get audited -> V(betat2); if audited -> (1-d_steal)*V(betat2) + d_steal*V(betat3)
  !Case 3: iterm = 2 & ipast = 3 (caught stealing): Construct value function only with beta3 -no matter their behavior, they will have the corrupt mayor betas. The difference
  !        is that if they are not audited, netW = newsav*R and if they are audited then netW = newsav*R - fine_fun(stealing)

  !Calculating the emaxp [future value] 
  !the value function depends on log of public consumption acounting for the degree of rivalry: qconsvec = log(Q/P**eta)
  ln_relqcons = LOG(qcons_w_riv)
  new_fr_stolen = newsteal/fundsl
  
  !Find the point on the saving domain where we need to compute the approximated value function
  DO isav = 1, Nwealth
     IF (newsav <  wealth(isav)) THEN
        ub_sav = isav
        EXIT
     END IF
     ub_sav = Nwealth
  END DO

  DO ifine = 1, Nfines
     frac_fine_g = frac_fine(ifine)
     netW(ifine) = newsav - fine_fun(newsteal)
     DO isav = 1, Nwealth
        IF (netW(ifine) <  wealth(isav)) THEN
           ub_netW(ifine) = isav
           EXIT
        END IF
        ub_netW(ifine) = Nwealth
     END DO
  END DO
  
  IF (Npast > 1) THEN     
     betat1_t = betat1(:)
     betat2_t = betat2(:)
  ELSE
     betat1_t = betat1(:)
     betat2_t = betat1_t
  END IF

  !Constructing emaxp1: value function if future audit = 0
  !In this case, savings = newsav independently of stealing
  IF(MMg == MMability_steal_q) THEN
     
     !COMPUTE THE LINEAR APPROXIMATION
     IF (ub_sav == 1) THEN
        
        multi_slope = DBLE(CEILING(ABS((newsav - wealth(1)))/inc_wealth))*0.0001d0
        V_slope = (1d0+multi_slope)*(betat1_t(2) - betat1_t(1))/(wealth(2) - wealth(1))
        
     ELSE

        V_slope = (betat1_t(ub_sav) - betat1_t(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))

     END IF

     ! Since they will not be audited, we don't include newsteal
     emaxp1 = betat1_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+1) + betat1_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+2)*ability_g + &
          betat1_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+3)*ability_g**2d0 + betat1_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+5)*ln_relqcons + &
          V_slope*(newsav - wealth(ub_sav)) + betat1_t(ub_sav)

!!$     if (myrank == 0 .AND. flag_simulation == 1) write(*,'(a,1i4,7f16.8)') 'emaxp1', ub_sav, V_slope*(newsav - wealth(ub_sav)) + betat1_t(ub_sav), V_slope*(newsav - wealth(ub_sav)), betat1_t(ub_sav), V_slope, (newsav - wealth(ub_sav)), newsav, wealth(ub_sav)
        
  ELSE
     
     !COMPUTE THE LINEAR APPROXIMATION
     IF (ub_sav == 1) THEN
        
        multi_slope = DBLE(CEILING(ABS((newsav - wealth(1)))/inc_wealth))*0.0001d0
        V_slope = (1d0+multi_slope)*(betat1_t(2) - betat1_t(1))/(wealth(2) - wealth(1))
        
     ELSE
        
        V_slope = (betat1_t(ub_sav) - betat1_t(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
        
     END IF
     
     emaxp1 = V_slope*(newsav - wealth(ub_sav)) + betat1_t(ub_sav)

  END IF
     
  !Constructing emaxp2: value function if future audit=1; in this case savings = netW
  IF(MMg == MMability_steal_q) THEN

     DO ifine = 1, Nfines
 
        !COMPUTE THE LINEAR APPROXIMATION
        IF (ub_netW(ifine) == 1) THEN
           
           multi_slope = DBLE(CEILING(ABS((netW(ifine) - wealth(1)))/inc_wealth))*0.0001d0
           V_slope = (1d0+multi_slope)*(betat2_t(2) - betat2_t(1))/(wealth(2) - wealth(1))
           
        ELSE
           
           V_slope = (betat2_t(ub_netW(ifine)) - betat2_t(ub_netW(ifine)-1))/(wealth(ub_netW(ifine)) - wealth(ub_netW(ifine)-1))
           
        END IF
        
        ! Since they will be audited we include newsteal without any dummy for being audited
        emaxp2(ifine) = betat2_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+1) + betat2_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+2)*ability_g + &
             betat2_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+3)*ability_g**2d0 + betat2_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+4)*new_fr_stolen + &
             betat2_t(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+5)*ln_relqcons + V_slope*(netW(ifine) - wealth(ub_netW(ifine))) + betat2_t(ub_netW(ifine))

        emaxp2(ifine) = emaxp2(ifine)*pfine(ifine)
        
     END DO

  ELSE

     DO ifine = 1, Nfines
        
        !COMPUTE THE LINEAR APPROXIMATION
        IF (ub_netW(ifine) == 1) THEN
           
           multi_slope = DBLE(CEILING(ABS((netW(ifine) - wealth(1)))/inc_wealth))*0.0001d0
           V_slope = (1d0+multi_slope)*(betat2_t(2) - betat2_t(1))/(wealth(2) - wealth(1))
           
        ELSE
           
           V_slope = (betat2_t(ub_netW(ifine)) - betat2_t(ub_netW(ifine)-1))/(wealth(ub_netW(ifine)) - wealth(ub_netW(ifine)-1))
           
        END IF
        
        emaxp2(ifine) = V_slope*(netW(ifine) - wealth(ub_netW(ifine))) + betat2_t(ub_netW(ifine))

        emaxp2(ifine) = emaxp2(ifine)*pfine(ifine)
        
     END DO

  END IF

  !Constructing total future value: pr(au=0)*f(betat1_t())+ pr(au=1)*[(1-dum_steal)*f(betat2_t()) + dum_steal*f(betat3_t())]
  emaxp = 0d0
  emaxp = paudit_g*emaxp1 + (1d0-paudit_g)*SUM(emaxp2(:))

  ! If flag_no_run == 1 .AND. new_fr_stolen > 0d0 the incumbent cannot run
  IF ((flag_simulation == 0 .AND. flag_no_run == 1 .AND. new_fr_stolen > 0d0 .AND. iterml < Nnterms) &
       .OR. (flag_simulation == 1 .AND. flag_no_run == 1 .AND. new_fr_stolen > 0d0 .AND. iterml < Nnterms .AND. timel < Retage -1)) THEN

     betat3_t = betat3(:)

     DO ifine = 1, Nfines
        
        !COMPUTE THE LINEAR APPROXIMATION
        IF (ub_netW(ifine) == 1) THEN
           
           multi_slope = DBLE(CEILING(ABS((netW(ifine) - wealth(1)))/inc_wealth))*0.0001d0
           V_slope = (1d0+multi_slope)*(betat3_t(2) - betat3_t(1))/(wealth(2) - wealth(1))
           
        ELSE
           
           V_slope = (betat3_t(ub_netW(ifine)) - betat3_t(ub_netW(ifine)-1))/(wealth(ub_netW(ifine)) - wealth(ub_netW(ifine)-1))
           
        END IF
        
        emaxp3(ifine) = V_slope*(netW(ifine) - wealth(ub_netW(ifine))) + betat3_t(ub_netW(ifine))

        emaxp3(ifine) = emaxp3(ifine)*pfine(ifine)
        
     END DO

     !Constructing total future value: pr(au=0)*f(betat1_t())+ pr(au=1)*[(1-dum_steal)*f(betat2_t()) + dum_steal*f(betat3_t())]
     emaxp = 0d0
     emaxp = paudit_g*emaxp1 + (1d0-paudit_g)*SUM(emaxp3(:))

  END IF

  ! If flag_no_run == 1 .AND. new_fr_stolen > 0d0 the incumbent cannot run
  IF ((flag_simulation == 0 .AND. flag_norun_prob_policy == 1 .AND. new_fr_stolen > l_fr_stolen_x .AND. iterml < Nnterms) &
       .OR. (flag_simulation == 1 .AND. flag_norun_prob_policy == 1 .AND. new_fr_stolen > l_fr_stolen_x .AND. iterml < Nnterms .AND. timel < Retage -1)) THEN

     betat3_t = betat3(:)

     DO ifine = 1, Nfines
        
        !COMPUTE THE LINEAR APPROXIMATION
        IF (ub_netW(ifine) == 1) THEN
           
           multi_slope = DBLE(CEILING(ABS((netW(ifine) - wealth(1)))/inc_wealth))*0.0001d0
           V_slope = (1d0+multi_slope)*(betat3_t(2) - betat3_t(1))/(wealth(2) - wealth(1))
           
        ELSE
           
           V_slope = (betat3_t(ub_netW(ifine)) - betat3_t(ub_netW(ifine)-1))/(wealth(ub_netW(ifine)) - wealth(ub_netW(ifine)-1))
           
        END IF
        
        emaxp3(ifine) = V_slope*(netW(ifine) - wealth(ub_netW(ifine))) + betat3_t(ub_netW(ifine))

        emaxp3(ifine) = emaxp3(ifine)*pfine(ifine)
        
     END DO

     !Constructing total future value: pr(au=0)*f(betat1_t())+ pr(au=1)*[(1-dum_steal)*f(betat2_t()) + dum_steal*f(betat3_t())]
     emaxp = 0d0
     emaxp = paudit_g*emaxp1 + (1d0-paudit_g)*(SUM(emaxp2(:))*(1d0 - prob_conv_norun)+SUM(emaxp3(:))*prob_conv_norun)

  END IF


  
  !Calculating total value [utility+future value]
  vemax = util + betadf*emaxp

  IF(vemax > vemaxoptg) THEN
     vemaxoptg = vemax
     pconsoptg = pcons
     !qconsoptg is public consumption in levels, it DOES NOT takes into account the degree of rivalry
     qconsoptg = qcons
     newsavoptg = newsav
     stealoptg = newsteal
     utiloptg = util
     emaxoptg = emaxp
     zzpuoptg = zzpu
  END IF

  fval = -vemax

!!$  if (myrank == 0 .AND. flag_simulation == 1) THEN
!!$  if (myrank == 0 .AND. flag_here == 124) then
!!$  if (flag_here == 246) then
!!$     write(*,'(a,3i4,13f15.8)') 'mayor sol 0', MMg, iterml, Nnterms, vemax, vemaxoptg, newsteal, newsav, pcons, qcons_w_riv, util, zzpu, emaxp, emaxp1, emaxp2, emaxp3, paudit_g
!!$     write(*,'(a,16f14.8)') 'mayor sol 1', betat1_t 
!!$     write(*,'(a,16f14.8)') 'mayor sol 2', betat2_t 
!!$     write(*,'(a,16f14.8)') 'mayor sol 3', betat3_t 
!!$     write(*,*) ''
!!$  end if
  
END SUBROUTINE mayorsolution
